Modulational Instability in Nonlinearity-Managed Optical Media 
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We investigate analytically, numerically, and experimentally the modulational instability in a 
layered, cubically-nonlinear (Kerr) optical medium that consists of alternating layers of glass and 
air. We model this setting using a nonlinear Schrodinger (NLS) equation with a piecewise constant 
nonlinearity coefficient and conduct a theoretical analysis of its linear stability, obtaining a Kronig- 
Penney equation whose forbidden bands correspond to the modulationally unstable regimes. We 
find very good quantitative agreement between the theoretical analysis of the Kronig-Penney equa- 
tion, numerical simulations of the NLS equation, and the experimental results for the modulational 
instability. Because of the periodicity in the evolution variable arising from the layered medium, we 
find multiple instability regions rather than just the one that would occur in uniform media. 

PACS numbers: 05.45.Yv, 42.65.Sf, 42.65.Tg, 42.65.-k 
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I. INTRODUCTION 

The modulational instability (MI) is a destabilization 
mechanism for plane waves that results from the inter- 
play between nonlinear and dispersive effects [H, [!]■ It 
leads to delocalization in momentum space and, in turn, 
to localization in position space and the formation of lo- 
calized (solitary- wave) structures. The MI arises in many 
physical settings, including fluid dynamics (where it is 
also called the "Benjaniin-Feir instability") nonlin- 
ear optics [1,111, [1], plasma physics and other areas. 
Recently, it has also taken center stage in atomic physics 
in studies of atomic Bose-Einstein condensates (BECs) 
0,113, [HE [II (see also [H for a review). In all of the 
above settings, the MI is one of the principal mechanisms 
leading to the emergence of localized "coherent nonlin- 
ear structures" and the formation of bright solitary waves 
(and trains of solitary waves). 

The MI was originally analyzed in uniform me- 
dia, predominantly in the framework of the nonlinear 
Schrodinger equation (NLS) [IB], where a focusing non- 
linearity leads to MI for sufficiently large plane- wave am- 
plitudes (for a given wavenumber) or sufficiently small 
wavenumbers (for a given amplitude) ^3] . More recently, 
several very interesting experimentally relevant settings 
with (temporally and/or spatially) nonuniform media 
have emerged. Research in this direction includes the 
experimental observation of bright matter-wave soliton 
trains in BECs [l6l |. which were induced by a tempo- 
ral change of the interatomic interaction from repulsive 
to attractive through Feshbach resonances. Subsequent 
theoretical work demonstrated how this effective change 
of the nonlinearity from defocusing to focusing leads 
to the onset of MI and the formation of soliton trains 
[Til . [T^ . [ist . Such soliton trains can also be induced in 
optical settings, as has been demonstrated, for example, 
in the context of birefringent dispersive media [itI . [18| . 



The theme of the present article is the analysis of the 
MI in a setting with periodic nonuniformities. In the con- 
text of nonlinear science, such studies date back to the 
Fermi-Pasta-Ulam problem [3, [I3| and have continued 
ubiquitously in numerous disciplines. Periodic nonuni- 
formities arise not only in the physical sciences (optics, 
atomic physics, solid-state physics, etc.) but also in biol- 
ogy in, for example, DNA double-strand dynamics [2Q|] . 
In optics, prominent areas in which periodicity takes cen- 
ter stage include the study of photonic crystals [IH and 
optically-induced lattices in photorefractive crystals [H] . 
Additionally, there have been recent experimental ob- 
servations of the MI in spatially periodic optical media 
(waveguide arrays) [13] ■ In solid-state physics, periodic- 
ity has been prevalent in the study of superconducting 
Josephson junctions [23 |. among other topics. Similar 
studies have been reported in atomic physics in, for ex- 
ample, BECs confined in spatially periodic optical po- 
tentials (so-called "optical lattices" ) [1, [13, [H, [ll] . A 
key feature in such settings is that MI can even occur for 
defocusing nonlinearities for certain wavenumber bands, 
as was originally suggested in Ref. 

While the aforementioned results pertain to spatially 
inhomogeneous settings, in which the periodicity lies in 
the transverse dimensions, we discuss in this article the 
theoretical analysis and the first experimental realization 
of MI in a setting that is periodic in the evolution vari- 
able. In the optical setting discussed here, this variable 
describes the propagation distance (it represents time in 
the framework of BECs). Such settings were initially 
proposed in the context of optical fiber communications, 
through so-called "dispersion management" techniques 
(which induce periodic changes in the group- velocity dis- 
persion) [H, [53, [13, [U, [13]), and have since also been 
studied for "nonlinearity management" (in which the 
Kerr nonlinearity is periodic in the propagation variable) 
[13, [H, [13, [1^- They have also been investigated in 
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BECs through so-called "Feshbach Resonance Manage- 
ment" (FRM) schemes [33, |38|, IM ES IM 113, El, 13 (see 
also the review p5l|). 

In the present paper, which provides a detailed de- 
scription of work we reported in a recent Letter (46j . we 
discuss an experimental investigation of MI in a layered 
optical medium and its quantitative comparison with an- 
alytical and numerical results. The medium consists of 
alternating layers of glass and air, through which the 
optical waves traverse. We model the dynamics of our 
experiment using an NLS equation with piecewise con- 
stant nonlinearity coefficients and a dissipation mecha- 
nism that accounts for reflective losses at the glass-air 
interfaces. We obtain very good quantitative agreement 
between our theoretical analysis, numerical simulations, 
and the experimental results. The tractability of the the- 
oretical analysis in our setting stems from the piecewise- 
constant nature of the material coefficients. In particular, 
this leads to an ordinary differential equation (ODE) in 
the linear stability analysis of the plane waves that is rem- 
iniscent of the Kronig-Penney (KP) model of solid-state 
physics [131 . This model is a special case of a Hill equa- 
tion [i^, whose coefficients are periodic in the evolution 
variable, leading to a band structure for its eigenvalues. 
In the present context, the forbidden bands of the KP 
model correspond to the instability regimes. This ob- 
servation allows us to compute the bands of modulation- 
ally unstable wavenumbers of perturbation (on top of the 
plane-wave solution) semi-analytically and consequently 
to compare our experimental findings not only with the 
numerical simulations of the NLS equation but also with 
the theoretical analysis of the KP model. Similar analy- 
ses have been performed, from a theoretical point of view, 
for dispersion-managed optical fibers [H, [23, [3l| . 

The remainder of this paper is organized as follows. 
First, we present our theoretical model governing pulse 
propagation in layered Kerr media and we analyze it 
mathematically. We then discuss our numerical proce- 
dures and experimental setup, present our main results, 
and provide additional discussion on several technical 
points. Finally, we present our conclusions and some 
future challenges. 



we derive the following NLS equation: 

i-^^--\7lu-\u\^u, 0<C<1 (glass). 



du 



1 n, 



(1) 



,(2) 



^^-k^i"-^!"!"' 1<C<L (air). 



(1) 



In Eq. ([T]), space is rescaled by the wavenumber k'^^^ — 
2nn'^^/X (where A is the source wavelength); that is, 
(C, '7, C) = fc^^"* X {^tDtz)- The complex field u is given 
by u = {n"2^ /rS^^Y^'^E, where E is the electric field en- 
velope. The superscript [j) denotes the medium, with 
j — 1 for glass and j — 2 for air. As is well-known, the 
linear parts of the refractive indices of glass and air 
take the values nl^^ = 1.5 and nQ^-* = 1, respectively. 
The nonlinear parts (i.e., the Kerr coefficients), are 



,(1) 



= 3.2 X 10-i« cmVW and nf^ = 3.2 x IQ-^^ cm^/W 

One can also incorporate the transmission losses at 
each slide into the system ([Ij. It is then written 



.du 



1 



D{C)V^u~NiC)\u\^u-ijiOu, (2) 



where D{C) = 1, N{C) = 1 in glass and £>(() n[,^V4^\ 
N{Q — in air. The last term in Eq. ([2|) de- 

scribes the transmission losses at each slide. The loss 
rate 7(C) is given by 



M 



y(C) = a5]<5(C-Cn) 



(3) 



where M is the number of glass-air interfaces at which 
losses occur and C„ is the location of the n-th interface. 
The prefactor a is determined by the constraint that the 
power P after an interface is a factor r (which for our 
experiments is typically 0.99) times the power before the 
interface. Because dP/dC, = — 27P, the parameter a sat- 
isfies the equation exp(— 2a) = r. 

To examine the onset of MI, we consider plane- wave 
solutions of Eq. ([1]), which are uniform in the tranverse 
spatial variables (^ and 77). Using the transformation 



II. THEORETICAL ANALYSIS 



A. Analytical results 



u exp 



we obtain the equation 



C 



7(C'X 



(4) 



Our theoretical model for the propagation of the opti- 
cal beam in the layered nonlinear medium incorporates 
the dominant dispersive and Kerr effects for each of the 
two media. Accordingly, by employing the slowly- varying 
envelope approximation to the Maxwell equations |49|], 



dC 



liC)dC 



\v\ V, 



(5) 



where we note that the Laplacian term in Eq. ([T]) is 
identically zero. Transforming to polar coordinates, v = 



3 



Re^^ , we subsequently obtain 
.dR de 



dR 



-iV(C)exp 



M 



dC 



= 0, 



(6) 



where i?(C') is the Heaviside step function. Equation 
impUes R — R{0) = Aq and 



d0 
dC 



M 



-2a^iJ(C-Cn) 



This yields an analytical expression for the plane-wave 
solutions of Eq. 0, 



(8) 



where r(C') = exp(-2 /'■' 7(C)rfC)- 

To perform a linear stability analysis of these plane 
waves, we consider a spatial perturbation of ([8]) given by 



u = uo(C) [1 + w(C) cos(fc^i^) cos(fc^77)] 



(9) 



where w is a small (Fourier-mode) perturbation with 
wavevector {k^,k^). We use the notation w = ew and 
insert the expression ^ into Eq. ^ to obtain 

«^(C)[1 + ew(C) cos(fc5^) cos(fc^77)] 

+ ieuoiO^iC) cos(%0 cos(fc,,?7) + 0{e^) 

= -^sMODiOhkl - k'^)w{C) cos(%C) cos(fc^?7) 
- *7(C)wo(C)[l + u){C) cos(fc^^) coa{kjji])] 

- iV(C)|wo{C)l'wo(C) [1 + 2ew{0 cos(fc^e) cos(fc„7,)] x 

X [1 + ew*{0 cos(%0 cos(fc^r;)] + 0{s^) , (10) 

where w* denotes the complex conjugate of li. We then 
equate the terms in (jlOp order by order in powers of e. 
By construction, the 0(1) terms cancel out. The 0{e) 
terms give 

'^iCMO + ^MO^iO = liki + k^,)DiOuo{OwiO 

- z7(C)«o(C)t&(C) - NiO\uo{(:)\^uo{0{2w + w*) . 



(11) 



Dividing both sides of (fTT|) by uq{() and using the equa- 
tion 



. 1 duo 

Wo C dC 



-iV(C)|uo(C)p-n(C), 



yields 



dw , , 1 

- N{0\uo{0?w{0 - N{0\uo{0?w*{0 . (12) 



^^iC) = ^ikl + kl)D{C)w{0 



Decomposing w (and hence w) into real and imaginary 
parts, w = F + iB, we obtain a linear system of ODEs, 



dC 2 ' 



(13) 



dB 



1 



^ = -2 P^(C) + 2A^(C)K(C)l']^, 



where k"^ = fc| -|- fc^. We rewrite the above system as a 
single second-order equation. 



(7) 'EL 



d^F 
dC 



+ 



1 dD dF 
\~k'D{Cr 



n{C)Pd{C)\uo{c)\' 



F.(14) 



The transformation F = gD^/^ then yields a Hill equa- 
tion, 



dC^ 



s{0 



1 



dD 



ADicr \dc 



(C) 



1 1 d^D 

2U(x~dC^ 



(C) 



s{0 - 2p^(C) 



27V(C)K(C)P-2^'^^^) 



5, 
(15) 



Hereafter, we make some simplifying assumptions that 
we will justify based both on the experimental system pa- 
rameters and on a detailed comparison of our results with 
the numerical simulations and the experimental findings. 
In particular, in investigating Eq. ([T5)) . we ignore the 
losses at the interfaces so that its solution is periodic 
and Bloch's theorem can be employed. Given that the 
transmission for each slide is approximately 99%, this is a 
very good approximation. Additionally, while it is possi- 
ble to analyze Eq. directly (if losses are ignored), 
it is sufRcient in modeling our experiments to exploit 
the weak variation of D{Q and substitute D{() with its 
average. In this case, the transformation to Eq. (fT5|) 
is no longer necessary, and we proceed using Eq. (jl4p . 
As will be discussed below, this assumption has hardly 
any effect on the results from Eq. p^ . Under this ad- 
ditional simplification, Eq. is a Hill equation 
equivalent to the well-known Kronig-Penney model from 
solid-state physics (originally aimed as a prototypical de- 
scription of an electron moving through a crystal lattice 
in a solid) [47|] for the piecewise-constant nonlinearity co- 
efficient under consideration. As discussed in Ref. [isj . 
such a periodic potential allows one to use Bloch's the- 
orem to obtain an analytical solution in both glass and 
air. Consequently, one may write [47l l5l| 



F(C + £)=exp -jc^Z F(C), 



(16) 



where to is the Floquet multipler. The analytical solution 
for F{C) is 

F = aie'-'"^ + bie-''^^ , < C < ^ (glass) , 
F = aae'''^'^ + b2e~''^^ , 1<C<L (air) , 
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where sf = k'^ D^'^k"" D<^^^ / i - TV^iVoP), sj = 
/t2D(2)(fc2£,(2)/4_^(2)|y^|2)^ g^^^ ^ ^l^g Floquet mul- 
tiplier. Because we are ignoring losses at the interfaces, 
|uoP = l^oP in the expressions for si and S2. 
The continuity of F and dF/dC, at the glass-air bound- 



aries (i.e., at C = ' and ( = L) leads to matching con- 
ditions that are used to determine the constants of in- 
tegration ai, 02, &i, and 62- This yields the following 
homogeneous 4x4 system of equations: 



/ 1 



1 



e 

-Sie 



-1 

-S2 



-1 
S2 



_g-ja;Lg- 



526 



-iS2 (i- 
-iS2{l- 




0, 



r 



which possesses nontrivial solutions if and only if the de- 
terminant of the matrix vanishes. This gives the follow- 
ing solvability condition for lo: 

g2 g2 ^ _ 

cos(cj_L) = sin(si^) sin[s2(i — I)] 

2siS2 

+ cos(si[) cos[s2(i - 0] = G(k) , (17) 

similar to that obtained for plane-wave solutions of the 
NLS equation with piecewise constant dispersion man- 
agement [15. Because of the functional form of the 
left-hand-side of Eq. P7|) . real solutions for the Flo- 
quet exponent lo exist if and only if 10(^)1 < 1. Given 
the form of the solution for the perturbation F{(), this 
case corresponds to modulationally stable wavenumbers 
k that exhibit oscillatory behavior. On the other hand, 
for |G'(fc)| > 1, the solutions of Eq. (fTT]) are imaginary, 
leading to an exponential growth in the real part of the 
perturbation F(^) which, in turn, indicates that such 
wavenumbers are modulationally unstable. Hence, the 
analogy to the original Kronig-Penney problem shows 
that the allowable energy zones are the ones correspond- 
ing to stable wavenumbers, whereas the forbidden energy 
zones are the ones associated with MI. 



B. Numerical setup 

We perform direct numerical simulations of Eq. ^ 
using experimentally-determined parameters. To con- 
firm the validity of our simulations, we use two different 
algorithms: a beam-propagation, split-step code and a 
code using finite differences in the spatial variables and a 
fourth-order Runge-Kutta algorithm in the propagation 
direction. We used the latter method to generate all of 
the figures shown in this paper and have verified that we 
obtain the same results using the former algorithm. 

For the depicted figures, we used a grid with 1000 
spatial nodes and scaled the spatial step size by the 
wavenumber for each simulation to ensure periodic 
boundary conditions in a domain encompassing 15 pe- 
riods of the periodic initial condition [see, in particular, 



Eq. We obtained the same results upon varying the 
number of grid points (up to 4000) and the number of 
periods. The typical grid spacing and evolution-variable 
step were approximately 0.3/xm and 0.04^m, respectively. 

In conducting our numerical experiments, we make two 
additional simplifications. First, motivated by the exper- 
iment, we consider in our numerical simulations the one- 
dimensional dynamics along the direction of the modu- 
lation. That is, we use A:^ = and vary k^. Accord- 
ingly, we convert the experimental two-dimensional in- 
terference patterns recorded on the CCD camera (see 
the next section for the specifications) to one-dimensional 
ones by integrating along the direction orthogonal to the 
modulation. Second, we assume that the modulational 
dynamics of the (weakly decaying) central part of the 
Gaussian beam of the experiment is similar to that of 
a plane wave with the same intensity. We tested both 
of these assumptions and confirmed them both a priori 
through the dynamical evolution of our experimental and 
numerical results and a posteriori through their quantita- 
tive comparison. Consequently, the initial wavefunctions 
in the numerical experiments take the form 

u = Ao + eoe''^« , (18) 

where ^0 = 's/Z); and /q denotes the intensities used 
experimentally (see the discussion below). We also use 
the experimental perturbation value of eq = Aq/10. We 
follow (and report) the subsequent evolution in both real 
and Fourier space. 

III. EXPERIMENTS VS. THEORY 

A. Experimental setup 

In our experiments, which are shown schematically in 
Fig. [Tl we use an amplified Titanium:Sapphire laser to 
generate 150- femtosecond pulses with an energy of 2 mJ 
at a wavelength of A = 800 nm. The beam profile is ap- 
proximately Gaussian with a full-width at half-maximum 
(FWHM) of 1.5 mm. The laser pulses are split into a 
pump and a reference using a beam splitter (BSl), with 
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FIG. 1; Experimental setup. BSl and BS2 are beam splitters, 
DL is a variable delay line, Ml and M2 are mirrors, NLM is 
the layered nonlinear medium, and LI and L2 are lenses. 



most of the energy in the pump pulse. After synchro- 
nization with a variable delay line (DL), the two pulses 
are recombined at a second beam splitter (BS2) and sent 
to the periodic nonlinear medium (NLM). This configu- 
ration allows us to control the relative angle in the prop- 
agation of the two beams. 

The reference introduces a sinusoidal modulation in 
the intensity (i.e., an interference pattern), with the pe- 
riod determined by the relative angle between the two 
beams. We carefully tune the angle of the reference by ro- 
tating BS2 so that the two beams overlap while propagat- 
ing through the NLM (at adjustable angles). The NLM 
consists of six 1 mm thick quartz microscope slides sepa- 
rated by air gaps. The glass slides have an anti-reflection 
coating to minimize the loss (the reflection from each in- 
terface is 1%). The loss due to back-reflections from the 
slides is included in our numerical simulations, and the 
effect of double reflections is negligible. In our experi- 
ments, we used structures with air gaps of 2.1 mm and 
3.1 mm. We image the intensity pattern after the NLM 
(at the output face of the last quartz slide) on a CCD 
camera (Pulnix TM-7EX) using two lenses (LI and L2) 
in a 4-F configuration with a magnification of A/ = 8. 
An image of the pump beam is used as a background 
and subtracted from the interference pattern to remove 
spatial nonuniformities that are not due to MI. The CCD 
camera captures the central region of the beam (0.6 mm 
X 0.8 mm). (Because the decay in this region of the 
Gaussian beam is weak, we approximate it as a plane 
wave in the theory and computations.) 

We record the intensity pattern at the output of the 
NLM for three different values of the pump intensity: 
Ipo = 9.0 X 10^ W/cm^, Jpi = 9.0 x 10^° W/cm^ and 
Ip2 = 1.3 X 10^^ W/cm^. In all three cases, the intensity 
of the reference beam is 1% of that of the pump. We 
measure the effect of the nonlinearity by comparing the 
output for high {Ipi and Ip2) versus low intensity (Ipo). 
For low intensity, the propagation is essentially linear. In 
the nonlinear regime, if the spatial frequency of the mod- 
ulation lies within an instability window, the amplitude 
of the reference wave will increase at the expense of the 
pump. 



B. Results 

The input field is given by Eq. p8|) . where Aq and 
eo are the amplitudes of the pump and reference beams, 
respectively, and |eoP <C |^oP- For hnear propagation 
(low pump intensity, Ipo), the intensity pattern at the 
output of the NLM is approximately the same as that at 
the input. That is, it is about 



loiO = |^o|' + ko|' + 2AoeoCos(%a. 



(19) 



For the nonlinear case (high pump intensity, Ipi and 
/p2), the amplitude of the waves changes and higher spa- 
tial harmonics are generated. The intensity at the output 
of the NLM is approximately 

iNLiO = |^i|^ + 2AieiCOs(fc5^) + 2Aie2Cos(2/ceO + --- , 

(20) 

where Ai and e„ {n — 1,2,---) are, respectively, the 
amplitudes of the pump beam and the n-th harmonic at 
the output of the NLM. The Fourier transform (FT) of 
Eq. jlQl) is 



FTilNL)^\A,\'S{f^) 
+ Aiei 

+ Ai€2 



5\k- 



27r 



^-1 

TT 



27r 

k 

TT 



(21) 



The FT peak height ratios (first-order:zeroth-order), 
''0 = ^o/Aq and ri = ei/Ai, are approximately equal 
to the amplitude ratios of the reference and pump waves. 
(For the experimental value of eg — Aq/1Q, the error 
introduced by this approximation is roughly 1%.). In 
the linear case, the value of tq remains constant. How- 
ever, in the nonlinear case, the value of ri depends on 
the length of the nonlinear medium. We thus use the 
ratio R = ri/rg as a diagnostic measure for both our 
experimental and numerical results, so that i? > 1 indi- 
cates growth of the perturbation. This measurement is 
equivalent to the ratio ri{( = ()/ri{( — 0) (where ( is 
the scaled length of the nonlinear medium), which com- 
pares the amplitude of the reference wave at the output 
to that at the input, but it is more robust experimen- 
tally because it accounts for other linear effects (such as 
the limited coherence length and spatial overlap of the 
pulses) that can affect the strength of the peaks in the 
Fourier transform. Therefore, the value of R reflects only 
the changes that are caused by the nonlinearity. In the 
numerical simulations, the peaks in the Fourier trans- 
form are sharp (one pixel), whereas they are broader in 
the experiments. Accordingly, when computing R from 
the experimental data, we use the area under the peaks 
instead of the peak value. 

Figure shows the ratio R{k) for the structure with 1 
mm slabs of glass sandwiching 2.1 mm air spacings, where 
k = k^riQ^^ is the sine of the angle between the pump and 
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FIG. 2: (Color online) Comparison of experimental (top), nu- 
merical (middle), and analytical (bottom) results for the 1 
mm glass-2.1 mm air configuration as a function of the di- 
mensionless wavenumber k. For the diagnostics R and \G\ 
(defined in the text), values larger than 1 correspond to MI. 

reference beams. There are two instability bands, quan- 
tified experimentally by i? > 1, within the measurement 
range. Similar to what has been shown computationally 
for dispersion- managed media [1^, [3l[ , the periodicity in 
the evolution variable from the layered ( "nonlinearity- 
managed") medium induces a second instability band. 
Note that only a single band occurs in uniform me- 
dia. The maximum growth of the perturbation in the 
first and second bands appear dX k — 6.0 x lO^'^ and 
k = 1.70 X 10^2^ with values of i? = 2.05 and R = 1.19, 
respectively. The increase in the modulation is clearly 
visible in the ID intensity patterns (see Fig. [3]) . As in- 
dicated by the middle and bottom panels of Fig. [2l the 
positions of the instability bands are in very good agree- 
ment with both numerical and theoretical (indicated by 
the forbidden zones with \G\ > 1) predictions. The nu- 
merical simulations typically show a stronger instability 
than the experimental measurement; this results from 
the three-dimensional nature of the experiment that is 
not captured in the simulation. In the experiment, the 
spatial and temporal overlap of the two beams decreases 
with increasing fc, which leads to weakening of the higher- 
order peaks. Additionally, temporal dispersion leads to a 
reduction in the aggregate strength of the nonlinearity in 
the experiment. The small- amplitude ripples that appear 
in the numerical simulation within the stable region re- 
sult from the finite number of periods in the propagation 
distance rather than from actual instabilities. We discuss 
this issue in further detail below. Despite this difference, 
we stress that the simulations successfully achieve our 
primary goal of quantitatively capturing the locations of 
the instability windows. 

Figure [Slja) shows the normalized experimental one- 
dimensional intensity pattern at the output of the NLM 
for Ip2 (dashed curve) and /po (solid curve). The three 
panels correspond to the cases k = 4.2 x 10""^ (top), 
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FIG. 3: (Color online) Comparison of experimental (a) and 
numerical (b) one-dimensional intensity patterns at the out- 
put [z = 16.5 mm) of the nonlinear medium with 1 mm 
glass-2.1 mm air for k = 0.0042 [top panels of (a) and (b)], 
k — 0.0126 (middle), and k — 0.0170 (bottom), corresponding 
to the first instability band, the following stable region, and 
the second instability band, respectively. The dashed curves 
are for high intensity (/p2) and the solid ones are for low in- 
tensity (/po). To facilitate comparisons between curves with 
different initial intensities, we scale the curves in this figure 
(to "arbitrary units") using the condition that the mean of 
\u{x, z)f is 1. 



which lies in the first instability band; k = 1.26 x 10^^ 
(middle), in the stable region separating the two forbid- 
den zones; and k = 1.70 x 10"'^ (bottom), which lies in 
the second instability band. The comparison of high and 
low intensity clearly shows the effect of the nonlinearity 
on the propagation. When the modulation of the input 
wave lies within the instability band, the amplitude of 
the modulation for the high intensity wave increases due 
to MI. Wc have observed an increase in the amplitude of 
the modulation both in the first and second instability 
bands, whereas the modulation remains practically un- 
changed in the stable region. We show the corresponding 
numerical intensity patterns in Fig. [3I^b). As with the 
experimental results, the MI in the first and third pan- 
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FIG. 4: (Color online) Experimental Fourier spectra (top 
panels) and numerical Fourier spectra (bottom panels) at the 
end of propagation {z = 16.5 mm) of the layered structure 
with 6 glass slides (of thickness 1 mm), each pair of which 
sandwiches 2.1 mm of air. The left panels are for k — 0.0042 
(first instability band), the center ones are for k — 0.0126 
(stable region), and the right ones are for k = 0.0170 (second 
instability band). The dashed curves are for high intensity 
(/p2) and the solid ones are for low intensity (Ipo)- 



els causes the peaks in the intensity pattern to become 
higher and narrower, whereas this is not the case for the 
middle panel of the figure. We note in passing that a 
small increase in amplitude is discernible in the middle 
panels of Figs.[Ha,b). This is associated with the ripples 
mentioned above and will be discussed more extensively 
below. 

In Fig. m we show the Fourier transforms of the exper- 
imental (top panels) and numerical (bottom panels) in- 
tensity patterns at the output of our layered medium (i.e., 
at z = 16.5 mm). Observe the appearance of higher spa- 
tial harmonics of the initial modulation in the regions of 
instability. Such harmonics correspond to the narrowing 
of the peaks in the spatial interference pattern. The pan- 
els depict the wavenumbers k = 0.0042 (left), k = 0.0126 
(center), and k = 0.0170 (right). The first-order peaks 
(the ones closest to fc = 0) correspond to the modula- 



FIG. 5: (Color online) Same as Fig.[2]but for the 1 mm glass- 
3.1 mm air configuration. Observe the presence of a third MI 
band. 



tion of the input beam and are present for both low and 
high intensity. For high intensity, additional peaks ap- 
pear at the higher harmonics for unstable wavenumbers 
of the modulation (fc = 0.0042 and k = 0.0170). We 
also observed this harmonic generation in the numeri- 
cal simulation (bottom panels), in good agreement with 
the experiments. In contrast, such harmonics are absent 
in the experimental results for fc within the modulation- 
ally stable regions (i? < 1 in Fig. [21 see the top center 
panel of Fig. [4]), again in agreement with the theoretical 
prediction. As with the ripples observed in Fig. [5] and 
the weak intensity amplification of Fig. [31 the appear- 
ance of weaker harmonics in the numerical simulation for 
fc = 0.0126 arises from the use of finitely many propaga- 
tion periods in the numerical simulation; these result in 
a weak amplification even in modulationally stable cases. 
As we explain in detail below, incorporating additional 
propagation periods in the numerical evolution distin- 
guishes with increasing clarity the modulationally stable 
and unstable regions. 

Figure [51 shows the experimental (top), numerical 
(middle) and theoretical (bottom) instability windows for 
the new structure, in which the 1 mm glass slides are 
sandwiched between 3.1 mm air windows. The longer 
spatial period in the structure results in a smaller spac- 
ing between the instability bands in Fourier space. Once 
again, we obtain good quantitative agreement between 
experiment, numerics, and theory with respect to the lo- 
cations of the instability bands. In the experiments, the 
peaks of the first two bands are at fc = 5.3 x lO^'^ and 
fc = 1.46 X 10-'^, with values of i? = 2.13 and R = 1.35, 
respectively; a third band appears near fc = 2.05 x 10^^, 
with R = 1.03. We note that because of the larger 
air gaps, the instability windows shift towards lower 
wavenumbers. Observe, however, that there is a slight 
disparity in the location of the very weak (in the experi- 
ments) third band between the analytical, numerical, and 
experimental results. This may be attributable to the 
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FIG. 6: (Color online) Same as Fig. [3] but at the end of the 
propagation (z — 21.5 mm) for the 1 mm glass-3.1mm air 
configuration. The wavenumbers in the experimental plots 
are k = 0.0042 [top panels of (a) and (b)], k = 0.0146 (mid- 
dle), and k = 0.0205 (bottom). The first two wavenumbers 
are the same in the numerical plots, but we use k = 0.0208, 
the peak of the third numerical band, for the leist plot. These 
wavenumbers occur, respectively, in the first instability band, 
the second instability band, and the third instability band. 
As before, the dashed curves are for high intensity (/p2), and 
the solid ones are for low intensity (Ipo)- 



very weak growth rate of the instability in conjunction 
with our quasi- ID approximation versus the fuUy three- 
dimensional spatio-temporal nature of the experiment. 

Figure [6] shows the normahzed one-dimensional inten- 
sity pattern at the output of the NLM for Ip2 (dashed 
curve) and Ipo (solid curve). For both the experimental 
[Fig. [DJa)] and theoretical [Fig. [Hb)] results, the three 
panels are for k = 0.0044 (top), which lies in the first in- 
stability band; 0.0146 (middle), which lies in the second 
instabihty band; and k = 0.0205 (bottom) [k = 0.0208 
for the numerical results] , which lies in the third instabil- 
ity band. Once again, we clearly observe an increase in 
the modulation depth in the instability bands, whereas 
for the modulationally stable case, such an increase is 
absent. In fact, the experiment may even suggest a cor- 



responding decrease, which does not appear in the nu- 
merical simulations. We believe this decrease results from 
the non-ideal conditions of the experiments — particularly 
because the modulation of the input beam is not purely 
sinusoidal and the beam itself is not a plane wave. The 
third band shows only a weak instability. Experimen- 
tally, the larger angle between the two beams (leading 
to the high wavenumbers of this band) reduces the in- 
teraction between them because of the limited temporal 
and spatial overlap, thereby decreasing the strength of 
the nonlinearity. 

In Fig. [71 we show the Fourier transforms of the exper- 
imental (top panels) and numerical (bottom panels) in- 
tensity patterns at the output of the layered Kerr medium 
with the 3.1 mm air gaps. The appearance of higher spa- 
tial harmonics of the initial modulation in the regions 
of instability is again evident both experimentally and 
numerically, especially in the first two modulationally 
unstable bands. For the third band, the phenomenon 
is very weak, as discussed above. From left to right, 
the panels depict the results for wavenumbers in the first 
{k = 0.0044), second {k = 0.0146), and third (fc = 0.0205 
in the experiments and k — 0.0208 in the numerical sim- 
ulations) unstable bands. 



C. Discussion 

In this section, we discuss the evolution of u(x, z) and 
elaborate on several of the points mentioned in previ- 
ous sections. In particular, we discuss the evolution of 
the diagnostic R = R{z), previously reported only at the 
output of the layered medium (in order to compare with 
experiments). We also discuss how the different intensi- 
ties (and hence different effective nonlinearities) influence 
our experimental results, in terms of the diagnostic R as 
a function of wavenumber. For that same dependence 
(i? as a function of wavenumber fc), we consider the rip- 
ples previously discussed for the stable region and how 
their relative amplitudes compared to the peak heights 
in the instability regions vanish as the propagation dis- 
tance (that is, the number of propagation periods) in- 
creases. Finally, we validate the assumption (used in our 
theoretical analysis) of substituting the weak variation 
in the dispersion by its average by comparing the di- 
rect evolution results between the true dynamics and the 
average-dispersion ones. 

In Fig. [51 we show contour plots of the intensity 
|u(a;,z)|^ for our numerical simulations (for the layered 
medium with 1 mm glass slides sandwiching 2.1 mm air 
gaps) of low-intensity and high-intensity initial wave- 
functions at fc = 0.0042 (first instability band), k — 
0.0126 (stable region), and k — 0.0170 (second insta- 
bility band). For the stable region, the nonlinearity has 
only a small effect on the evolution of the interference 
pattern. For propagation in the instability bands, we ob- 
serve an increase in modulation depth in both cases but 
with marked differences in the evolution. For the lower 
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FIG. 7: (Color online) Same as Fig. |4] but at the end of thi 
propagation (z — 21.5 mm) for the 1 mm glass-3.1mm ai: 
configuration. As before, the experimental results are showi 
in (a) and the numerical ones are shown in (b) . The left panel; 
show the results for wavenumber k = 0.0044 (first instabiliti 
band), the middles panels are for k = 0.0146 (second band) 
and the right panels are for the third band {k — 0.0205 for th( 
experiments and k = 0.0208 for the numerical simulations) 
As before, the dashed curves are for high intensity (/p2) ant 
the solid ones are for low intensity (Ipo)- 



wavenumber, the modulation grows continuously and thi 
orientation of the pattern remains fixed. The most in- 
teresting behavior is observed in the second instability 
band, where the evolution of the pattern is markedly dif- 
ferent. The modulation depth increases and decreases 
periodically, with a net increase after each period. The 
orientation of the pattern shifts dramatically between the 
different materials, possibly suggesting a periodic change 
in the relative angle between the two beams. 

We now discuss in greater detail the dynamical depen- 



FIG. 8: Evolution of the intensity \u{x,z)f (in arbitrary 
units) in layered Kerr media consisting of 1 mm glass-2.1 
mm air layers. We simulate initial wavefunctions of both 
low intensity (left) and high intensity (right) plane waves 
with a small modulation of wavenumber k — 0.0042 (top), 
k = 0.0126 (middle), and k = 0.0170 (bottom). For each of 
the numerical simulations, we consider domains with 15 pe- 
riods. The ticks on the horizontal axes indicate the glass-air 
interfaces, and the labeled values indicate the left edges of the 
glass slides. 
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FIG. 9: (Color online) Evolution of the MI diagnostic R = 
R(z) as a function of the propagation distance (in /im) for the 
configuration with 1 mm glass slides sandwiching 2.1 mm of 
air. The wavenumbers are k — 0.0042 (top), k — 0.0126 (mid- 
dle), and k = 0.0170 (bottom). The ticks on the horizontal 
axis indicate the glass-air interfaces, and the values labeled 
on the axis indicate the left edges of the glass slides. 
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FIG. 10: (Color online) Experimental measurement of R ver- 
sus the dimensionless wavenumber k. The top panel corre- 
sponds to the 1 mm glass-2.1 mm air configuration with inten- 
sity 7p2, the middle panel corresponds to the same structure 
with lower intensity Ipi, and the bottom panel corresponds 
to the structure with 1 mm glass-3.1 mm air and intensity 
Ip2- Observe in this last panel the leftward shift of the MI 
bands and the presence of a third band. 



dence of the diagnostic R on the propagation distance z 
(see Fig. For the case of 2.1 mm air gaps, we show 
examples in the first instability band (fc — 0.0042), the 
second instabihty band {k — 0.0170), and in the region 
between the two instability bands (fc = 0.0126). For the 
second unstable region, observe that R{z) oscillates but 
with an increasing amplitude. Note additionally that the 
peak of the oscillation is not at the boundary between the 
two layers but rather near the half-period point. We ob- 
serve similar features for the configuration with 3.1 mm 
air gaps. 

Figure [10] shows the experimentally-measured instabil- 
ity windows [R — R{k)] for the 2.1 mm air-1 mm glass 
configuration with an initial beam of high intensity Ip2 
(top), one with lower intensity Ipi (middle), and for the 
3.1 mm air-1 mm glass configuration with high inten- 
sity (bottom). For the lower intensity case, the peaks 
are smaller because of the weaker nonlinearity, in conso- 
nance with the theoretical prediction. For the structure 
with larger air gaps, the instability bands shift towards 
lower k. Moreover, as seen by Eq. ([1]), one can obtain 
more instability windows by increasing the widths of the 
air gaps further. However, this is very difficult to achieve 
in experiments. The interaction length of the two beams 
is limited by their spatial and temporal overlap. In order 
to study structures with larger gaps, the beam diameter 
and pulse duration must be increased while keeping the 
same value of the intensity, which in our case was not 
possible due to the limited pulse energy. 
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FIG. 11: (Color online) Numerical calculation of the instabil- 
ity bands (for the structure with 2.1 mm air gaps) for prop- 
agation with 6, 11, and 21 layers of glass (top, middle, and 
bottom panels, respectively). The heights of the ripples in 
the stable regions remain constant, whereas the heights of 
the peaks due to MI grow exponentially with propagation 
distance. 



We now examine the diagnostic R{k) for our direct nu- 
merical simulations as a function of propagation distance 
(number of propagation periods). Our aim is to explain 
the ripples in the stable regions in the middle panels of 
Figs. [2] and m We stress that these ripples do not corre- 
spond to an actual instability. Instead, they result from 
the fact that the simulations propagate over a finite num- 
ber of periods. We observed earlier that no such ripples 
appear in the instability windows computed using Bloch 
theory, which implicitly assumes that the number of pe- 
riods in the propagation direction is infinite. We thus ex- 
pect the prominence of such ripples to decrease for direct 
numerical simulations with more layers of glass and air. 
As shown by the numerical instability peaks in Fig. [11] 
for propagation distances with 6 (top panel), 11 (mid- 
dle), and 21 (bottom) layers of glass (and 5, 10, and 20 
sandwiched air gaps), this is indeed the case. In these 
numerical experiments, we decreased the amplitude of 
the input reference wave (eg) from 0.1 to 10^'^ to pre- 
vent a saturation in the growth of the modulation for 
the longer propagation distances. The simulations show 
that the amplitudes of the ripples remain essentially un- 
changed, whereas the MI peaks grow exponentially with 
distance, as expected from the theory. The position and 
periodicity of the ripples changes with the different prop- 
agation distances, which can be explained by the manner 
in which the growth of the reference beam evolves with 
propagation distance. 

Finally, we revisit (as promised earlier) the assumption 
in our mathematical analysis of a uniform, mean-valued 
dispersion in our study of MI. To do this, we performed 
direct numerical simulations in which the coefficient of 



11 



5.5 I r 




0.005 0.01 0.015 0.02 

k 



FIG. 12: (Color online) Instability windows for 1 mm glass- 
2.1 mm air configuration for piecewise constant D{() (solid 
curve) and D(C) = D = constant (dashed curve). 

the Laplacian was uniform (assuming the values of the 
weighted averages used in the theory) rather than piece- 
wise constant. For the 1 mm glass-2.1 mm air configura- 
tion, the mean of ^(C) is D = (1.5 + 2.1)/3.1 « 1.16. 
For the 1 mm glass-3.1 mm air configuration, it is 
D = (1.5 + 3.1)/4.1 « 1.12. As shown for the former 
configuration in Fig. [T^] (we obtained similar results for 
the latter one), such changes result in almost no differ- 
ences in the location of the instability windows and only 
small differences in the sizes of the instability peaks. We 
therefore assert that this provides an excellent approxi- 
mation for determining the locations of the modulation- 
ally unstable wavenumber windows. The controllability 
and validity of our approximations can therefore be used 
to explain the very good quantitative agreement that we 
observe between our analytical, numerical, and experi- 
mental results, especially in light of the fact that there 
are no free/adjustable parameters in our model. 

IV. CONCLUSIONS 

In the present work, we provided an experimental real- 
ization of the modulational instability (MI) in a medium 
that is periodic in the evolution variable. We also de- 
scribed the location of the instability bands quantitatively 
by investigating the Hill equation obtained from a linear 
stability analysis of plane-wave solutions of a nonlinear 
Schrodinger equation with piecewise constant nonlincar- 
ity coefficients. In this case, the Hill equation becomes 
the Kronig-Penney model and thereby provides a direct 
association of the MI bands with the forbidden energy 
zones of that model. One of the unique features of the 
periodic medium in this respect (which can also be seen 



in the dispersion-managed nonlinear Schrodinger equa- 
tion [H,!!^]) is the opening of additional MI bands, such 
as the second and third bands discussed in detail in the 
present work. The precise location of the bands is chiefly 
determined by the details of the periodicity in the evo- 
lution variable (the thickness of the glass slides and the 
width of the air gaps). We note in passing that our the- 
oretical result can be further analyzed in the limit in 
which the glass slides are much wider than the air gaps 
(or vice versa), in which case one can infer, for exam- 
ple, that for the higher zones (indexed by n), the zone 
width is inversely proportional to n (s^ . However, we did 
not pursue this case in detail, as it was not experimen- 
tally tractable. We compared our mathematical analysis 
for the modulationally unstable bands to both numerical 
and experimental findings (using Fourier-space diagnos- 
tics to elucidate the instability of the latter) and found 
very good agreement. We also clearly observed higher 
spatial harmonics for modulationally unstable beams, re- 
vealing another characteristic trait of MI. 

Additionally, the efficacy of layered Kerr media (and 
other instances of nonlinearity management) for exam- 
ining interesting and important nonlinear dynamics goes 
far beyond the present work on MI. In a recent paper 
js^ . we used a similar setup (alternating layers of glass 
and air) to stabilize an optical pulse using nonlinearity 
management, showing that it can potentially provide a 
lossless self-guiding mechanism. Because both air and 
glass are focusing media, collapse or dispersion cannot 
be entirely prevented. In fact, in our setting, the pres- 
ence of very weak dissipation always appears to favor the 
scenario of eventual dispersion of the pulse. Neverthe- 
less, this occurs for propagation distances that are an 
order of magnitude larger than the typical ones of uni- 
form media (and, moreover, the setting can be improved 
considerably). We also captured our experimental re- 
sults qualitatively (and, when appropriate, also quanti- 
tatively) by a (2 -|- l)-dimensional nonlinear Schrodinger 
equation with a piecewise constant nonlinearity coeffi- 
cient (with losses incorporated at the glass-air bound- 
aries). The very good agreement between theoretical, 
numerical, and experimental results both here and in the 
aforementioned previous work [3^ suggests a variety of 
extensions not only in the present setting of layered Kerr 
media but also in nonlinearity-managed Bose-Einstein 
condensates (whose mean-field dynamics are also gov- 
erned by nonlinear Schrodinger equations). In partic- 
ular, it will be very interesting to realize nonlinearity 
management in three-dimensional Bose-Einstein conden- 
sates and examine its impact on the stability of coherent 
structures — including both solitary waves and more com- 
plex entities such as vortices bearing topological charge 
(see, e.g., 

These recent developments (especially in optics but 
also Bose-Einstein condensation) underscore the inter- 
est in further developing and expanding the mathemat- 
ical theory pertinent to such settings. In particular, it 
would be of interest to examine the well-posedness of 
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such temporally-modulated settings (either at the level 
of the full model or at the level of its averaged variants), 
an avenue of research that is only starting to develop 

MM- 

Acknowledgements. We acknowledge support from the 



DARPA Center for Optofluidic Integration (DP), the 
Caltech Information Science and Technology initiative 
(MC, MAP), the Alexander von Humboldt-Foundation 
(MC) and NSF-DMS-0204585, NSF-DMS-0505663, NSF- 
DMS-0619492, and NSF-CAREER (PGK). 



[1] M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys. 65, 
851 (1993). 

[2] A. Scott (ed.). Encyclopedia of Nonlinear Science, Lec- 
tures in Applied Mathematics (Routledge, Taylor & 
Francis Group, New York, NY, 2005). 

[3] T. B. Benjamin and J. E. Feir, J. Fluid Mech. 27, 417 

(1967) . 

[4] G. P. Agrawal, Nonlinear Fiber Optics (Academic Press, 
San Diego, CA, 1995). 

[5] A. Hasegawa and Y. Kodama, Sohtons in Optical Com- 
munication (Clarendon Press, Oxford, UK, 1995). 

[6] L. A. Ostrovskii, Sov. Phys. JETP 24, 797 (1967). 

[7] A. Hasegawa, Phys. Rev. Lett. 24, 1165 (1970). 

[8] T. Taniuti and H. Washimi, Phys. Rev. Lett. 21, 209 

(1968) . 

[9] B. Wu and Q. Niu, Phys. Rev. A 64, 061603(R) (2001). 
[10] V. V. Konotop and M. Salerno, Phys. Rev. A 65, 

021602(R) (2002). 
[11] G. Theocharis, Z. Rapti, P. G. Kevrekidis, D. J. 

Frantzeskakis, and V. V. Konotop, Phys. Rev. A 67, 

063610 (2003). 

[12] L. Salasnich, A. Parola, and L. Reatto, Phys. Rev. Lett. 

91, 080405 (2003). 
[13] L. D. Carr and J. Brand, Phys. Rev. A 70, 033607 (2004). 
[14] P. G. Kevrekidis and D. J. Frantzeskakis, Mod. Phys. 

Lett. B 18, 173 (2004). 
[15] C. Sulem and P.-L. Sulem, The Nonlinear Schrddinger 

Equation: Self-Focusing and Wave Collapse, no. 139 in 

Applied Mathematical Sciences (Springer- Verlag, New 

York, NY, 1999). 
[16] K. E. Strecker, G. B. Partridge, A. G. Truscott, and R. G. 

Hulet, Nature 417, 150 (2002). 
[17] E. Seve, G. MiUot, S. Wabnitz, T. Sylvestre, and H. Mail- 

lotte, J. Opt. Soc. Am. B 16, 1642 (1999). 
[18] S. Wabnitz, Phys. Rev. A 38, 2018 (1988). 
[19] D. K. Campbell, P. Rosenau, and G. Zaslavsky, Chaos 

15, 015101 (2005). 
[20] M. Peyrard, Nonlinearity 17, Rl (2004). 
[21] M. Soljacic and J. D. Joannopoulos, Nature Materials 3, 

211 (2004). 

[22] N. K. Efremidis, S. Sears, D. N. Christodoulides, 
J. W. Fleischer, and M. Segev, Phys. Rev. E 66, 046602 
(2002). 

[23] J. Meier, G. I. Stegeman, D. N. Christodoulides, Y. Sil- 
berberg, R. Morandotti, H. Yang, G. Salamo, M. Sorel, 
and J. S. Aitchison, Phys. Rev. Lett. 92, 163902 (2004). 

[24] J. Pfeiffer, M. Schuster, A. A. Abdumalikov Jr., and 
A. V. Ustinov, Phys. Rev. Lett. 96, 034103 (2006). 

[25] F. S. Cataliotti, L. Fallani, F. Ferlaino, C. Fort, P. Mad- 
daloni, and M. Inguscio, New J. Phys. 5, 71 (2003). 

[26] A. Smerzi, A. Trombettoni, P. G. Kevrekidis, and A. R. 
Bishop, Phys. Rev. Lett. 89, 170402 (2002). 

[27] Y. S. Kivshar and M. Peyrard, Phys. Rev. A 46, 3198 
(1992). 



[28] J. C. Bronski and J. N. Kutz, Opt. Lett. 21, 937 (1996). 
[29] A. Kumar, A. Labruyere, and P. Tchofo Dinda, Opt. 

Commun. 219, 221 (2003). 
[30] B. A. Malomed, Progress in Optics 43, 71 (2002). 
[31] N. J. Smith and N. J. Doran, Opt. Lett. 21, 570 (1996). 
[32] S. K. Turitsyn, E. G. Shapiro, S. B. Medvedev, M. P. 

Fedoruk, and V. K. Mezentsev, Comptes Rendus Phys. 

4, 145 (2003). 

[33] L. Berge, V. Mezentsev, J. Rasmussen, P. Crhstiansen, 
and Y. Gaididei, Opt. Lett. 25, 1037 (2000). 

[34] M. Centurion, M. A. Porter, P. G. Kevrekidis, and 
D. Psaltis, Phys. Rev. Lett. 97, 033903 (2006). 

[35] F. Ilday and F. Wise, Opt. Lett. 19, 470 (2002). 

[36] I. Towers and B. A. Malomed, J. Opt. Soc. Am. B 19, 
537 (2002). 

[37] F. K. AbduUaev, J. G. Caputo, R. A. Kraenkel, and B. A. 

Malomed, Phys. Rev. A 67, 013605 (2003). 
[38] F. K. AbduUaev, S. A. Darmanyan, A. Kobyakov, and 

F. Lederer, Phys. Lett. A 220, 213 (1996). 
[39] F. K. AbduUaev, A. M. Kamchatnov, V. V. Konotop, 

and V. A. Brazhnyi, Phys. Rev. Lett. 90, 230402 (2003). 
[40] J. J. Garci'a-RipoU, V. M. Perez-Garci'a, and P. Torres, 

Phys. Rev. Lett. 83, 1715 (1999). 
[41] P. G. Kevrekidis, G. Theocharis, D. J. Frantzeskakis, and 

B. A. Malomed, Phys. Rev. Lett. 90, 230401 (2003). 
[42] D. D. Montesinos, V. M. Perez-Garcia, and P. J. Torres, 

Physica D 191, 193 (2004). 
[43] Z. Rapti, G. Theocharis, P. G. Kevrekidis, D. J. 

Frantzeskakis, and B. A. Malomed, Phys. Scr. T107, 27 

(2004) . 

[44] H. Saito and M. Ueda, Phys. Rev. Lett. 90, 040403 
(2003). 

[45] B. A. Malomed, Soliton Management in Periodic Systems 

(Springer- Verlag, New York, NY, 2006). 
[46] M. Centurion, M. A. Porter, Y. Pu, P. G. Kevrekidis, 

D. J. Frantzeskakis, and D. Psaltis, Phys. Rev. Lett. 97, 

234101 (2006). 

[47] C. Kittel, Introduction to Solid State Physics (John Wi- 
ley & Sons, Inc., New York, NY, 1986). 

[48] W. Magnus and S. Winkler, Hill's Equation (Interscience 
Publishers, New York, NY, 1966). 

[49] Y. S. Kivshar and G. P. Agrawal, Optical Solitons: From 
Fibers to Photonic Crystals (Academic Press, San Diego, 
California, 2003). 

[50] R. W. Boyd, Nonlinear Optics (Academic Press, San 
Diego, CA, 2003). 

[51] R. H. Rand Lecture Notes on Nonlinear Vi- 
brations, a free online book available at 
http: //www. tarn. Cornell . edu/randdocs/nlvibe52 .pdf 

(2005) . 

[52] I. Gol'dman and V. Krivchenkov, Problems in Quantum 
Mechanics (Dover, New York, NY, 1993). 

[53] V. V. Konotop and P. Pacciani, Phys. Rev. Lett. 94, 
240405 (2005). 



13 



[54] P. G. Kcvrokidis, D. Pclinovsky, and A. Stefanov, J. 
Phys. A: Math. Gen. 39, 479 (2006). 



